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Abstract. Automatic and semi-automatic magnetic resonance angiog- 
raphy (MRA) segmentation techniques can potentially save radiologists 
large amounts of time required for manual segmentation and can facili- 
tate further data analysis. The proposed MRA segmentation method uses 
a mathematical modeling technique which is well-suited to the compli- 
cated curve-like structure of blood vessels. We define the segmentation 
task as an energy minimization over all 3D curves and use a level set 
method to search for a solution. Our approach is an extension of previ- 
ous level set segmentation techniques to higher co-dimension. 



1 Introduction 

The high-level goal of this research is to develop computer vision techniques 
for the segmentation of medical images. Automatic and semi-automatic vision 
techniques can potentially assist clinicians in this task, saving them much of the 
time required to manually segment large data sets. Specifically, we consider the 
segmentation of volumetric vasculature images, such as the magnetic resonance 
angiography (MRA) image pictured in Figure 1. 

As shown here, blood vessels appear in MRA images as bright curve-like 
patterns which may be noisy and have gaps. What is shown is a "maximum 
intensity projection" . The data is a stack of slices where most areas are dark, 
but vessels tend to be bright. This stack is collapsed into a single image for 
viewing by performing a projection through the stack that assigns to each pixel 
in the projection the brightest voxel over all slices. This image shows projections 
along three orthogonal axes. 

Thresholding is one possible approach to this segmentation problem and 
works adequately on the larger vessels. The problem arises in detecting the small 
vessels, and that is the objective of our work. Thresholding cannot be used for 
the small vessels for several reasons. The voxels may have an intensity that is 
a combination of the intensities of vessels and background if the vessel is only 
partially inside the voxel. This sampling artifact is called partial voluming. Other 
imaging conditions can cause some background areas to be as bright as other 
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vessel areas, complicating threshold selection. Finally, the images are often noisy, 
and methods using local contextual information can be more robust. 

Our method uses the fact that the underlying structures in the image are 
indeed 3D curves and evolves an initial curve into the curves in the data (the 
vessels). In particular, we explore techniques based on the concept of mean cur- 
vature flow, or curve-shortening flow, from the field of differential geometry. 
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Fig. 1. Maximum intensity projection of a phase-contrast MRA image of blood vessels 
in the brain. 



2 Curvature Evolution Methods 



Mean curvature evolution schemes for segmentation, implemented with level set 
methods, have become an important approach in computer vision [5,10,11]. This 
approach uses partial differential equations to control the evolution. An overview 
to the superset of techniques using related partial differential equations can be 
found in [4]. The fundamental concepts from mathematics from which mean 
curvature schemes derive were explored several years earlier when smooth closed 
curves in 2D were proven to shrink to a point under mean curvature motion [8, 
9]. Evans and Spruck and Chen, Giga, and Goto independently framed mean 
curvature flow of any hypersurface as a level set problem and proved existence, 
uniqueness, and stability of viscosity solutions [7,6]. For application to image 
segmentation, a vector field was induced on the embedding space, so that the 
evolution could be controlled by an image gradient field or other image data. The 
same results of existence, uniqueness, and stability of viscosity solutions were 
obtained for the modified evolution equations for the case of planar curves, and 
experiments on real-world images demonstrated the effectiveness of the approach 
[3,5]. 

Curves evolving in the plane became surfaces evolving in space, called min- 
imal surfaces [5]. Although the theorem on planar curves shrinking to a point 



Co-Dimension 2 Geodesic Active Contours for MRA Segmentation 3 

could not be extended to the case of surfaces evolving in 3D, the existence, 
uniqueness, and stability results of the level set formalism held analogously to 
the 2D case. Thus the method was feasible for evolving both curves in 2D and 
surfaces in 3D. Beyond elegant mathematics, spectacular results on real- world 
data sets established the method as an important segmentation tool in both do- 
mains. One fundamental limitation to these schemes has been that they describe 
only the flow of hypersurfaces, i.e., surfaces of co-dimension 1. 

Altschuler and Grayson studied the problem of curve-shortening flow for 
3D curves [1], and Ambrosio and Soner generalized the level set technique to 
arbitrary manifolds in arbitrary dimension. They provided the analogous results 
and extended their level set evolution equation to account for an additional 
vector field induced on the space [2]. 

We herein present the first implementation of geodesic active contours in 
3D, based on Ambrosio and Soner's work. Specifically, our system uses these 
techniques for automatic segmentation of blood vessels in MRA images. The 
dimension of the manifold is 1, and its co-dimension is 2. 



3 Mean Curvature Flow 

Intuitively, mean curvature flow refers to some curve evolving in time so that at 
each point, the velocity vector normal to the curve is equal to the mean curvature 
vector. This concept is normally defined for arbitrary generic surfaces, but only 
curves are necessary for this paper, so we have restricted the definition. More 
formally, let C(t), t > be a family of curves in 3? 2 or 3? 3 , N the normal for a 
given orientation. That is, C is a curve, and t represents the "time" parameter 
or the index into the family of curves, not position. The mean curvature flow 
equation is then given by the vector equation 

C t = kN (1) 

with given initial curve C(0) — Co, n the curvature of the curve, and Ct the 
time derivative of the curve. Note that since we consider only ID curves here, 
as opposed to evolving surfaces, the mean curvature is just the usual curvature 
of the curve. This motion is also called "curve-shortening flow" since it is the 
solution, obtained by Euler-Lagrange equations, to the problem of minimizing 
curve length: 

min / \C'(p)\dp 
where p is the spatial parameter of the curve. 



4 Level Set Method for Planar Curves 

We give the basic idea of the level set method [12] to evolve a planar curve C. 
Define a function u : 3? 2 — > 3? so that C is a level-set of u. We follow the conven- 
tion that C is, in particular, the zero level set of u, although this choice is not 
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necessary for the method. The function u is now an implicit representation of 
the curve C . The advantages of this representation are that it is intrinsic (inde- 
pendent of parameterization) and that it is topologically flexible since different 
topologies of C are represented by the constant topology of u. Let Co be the 
initial curve. 

It is shown in [7] and [6] that evolving C according to 

C t = (3N (2) 

with initial condition C(-, 0) — Co(-) for any function j3, is equivalent to evolving 
u according to 

u t = /3|V«| (3) 

with initial condition «(-, 0) = Mo(") an d uq(Cq) = 0. 



' u = J 
u = O 



Fig. 2. Level sets of an embedding function u, for a closed curve in 5ft 2 



This result is independent of the choice of function u [7, 6]. As customary in 
the literature, we choose uq to be the signed distance function to the curve C 
(Figure 2). 



5 Level Set Method for Curves in Higher Codimension 

The level set evolution equations that follow were proven in [2]. They enable us 
to evolve space curves, with evolution driven by both mean curvature and image 
information. In the following discussion, C is a curve in 3D. 

5.1 Mean Curvature Flow 

Let v : 3? 3 — > [0, oo) be an auxiliary function whose zero level set is identically C, 
that is smooth near C, and such that Vf is non-zero outside C. For a nonzero 
vector q£ S", define 

T 
P n =I ^ 



q 'q| 2 
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Fig. 3. Evolving curves under mean curvature flow. The first three images show a circle 
shrinking to a point, and the last two images show a helix shrinking to its axis. 



as the projector onto the plane normal to q. Further define A(Vv(x, t), V 2 v{x, £)) 
as the smaller nonzero eigenvalue of P\r v V 2 vP^/ v . The level set evolution equa- 
tion is then 

v t = \{Vv{x,t),V 2 v{x,t)). 

That is, this evolution is equivalent to evolving C according to C< = kN in the 
sense that C is the zero level set of v throughout the evolution. 

Figure 3 demonstrates this evolution. As discussed above, a circle shrinks to 
a point under mean curvature motion. Under this motion, a helix evolves into 
its axis. 



5.2 Incorporation of Vector Field 

This section discusses the situation where there is an underlying vector field driv- 
ing the evolution, in combination with the curvature term. Assume the desired 
evolution equation is of the form 

C t = nN - lid, 

where 77 is the projection operator onto the normal space of C (which is a vector 
space of dimension 2) and d is a given vector field in 3? 3 . The evolution equation 
for the embedding space then becomes 

v t = \(Vv,V 2 v) + Vied. 




Fig. 4. (a) The tangent to C at p, the normal plane, the image-based vector, and its 
projection onto the normal plane, (b) e-level set method. 
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Fig. 5. Evolving helix under mean curvature flow with additional vector field: target 
curve, initial level set, level set after evolution with endpoints constrained. 



5.3 3D Image Segmentation 

For the case of ID structures in 3D images, we wish to minimize 



g(\VI(C(p))\)\C'(p)\dp 



where C(p) : [0, 1] ->■ 3? 3 is the ID curve, / : [0, a] x [0, b] x [0, c] -*■ [0, oo) 
is the image, and g : [0, oo) — > 3? + is a strictly decreasing function such that 
g(r) — ► as r — > oo (analogous to [5]). For our current implementation, we use 
g(r) = exp(-r) because it works well in practice. Another common choice is 
<7(|V/|) = 1+ |v/i2 ■ By computing the Euler-Lagrange equations, we find that 
the curve evolution equation is 



q' VI 

c t = K N-Ln^ JW 



(4) 



where H is the Hessian of the intensity function. The second term in the above 
equation is illustrated in Figure 4(a). That is, 



d= 9 -H- 



VI 



VI 



so the equation for the embedding space is 



vt = \(Vv(x, t),V 2 v(x, t)) + ^Vv(x, t) ■ H— — 

9 I VI 



(5) 



Thus, Ambrosio and Soner's work has provided the basis for the use of mean 
curvature flow and level set methods to segment ID structures in 3D. Figure 5 
illustrates how underlying image information can attract the evolving tube. The 
underlying volumetric image data is shown, as a maximum intensity projection, 
in the first image. This volume was generated by drawing a cosine curve in 
the volume, then smoothing with a Gaussian filter. The second image shows the 
initial curve, a helix. The result of the evolution is shown in the rightmost image. 



6 MRA Segmentation System 

This section describes our system for segmentation of vessels from MRA using the 
described level set method. A flowchart is shown in Figure 6. We discuss issues 
that have arisen in converting the theory above to practice for this application. 
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Fig. 6. Overview of segmentation algorithm. 



e-Level Set Method: Since the projection operator P q is defined only for 
non-zero vectors q, the method is undefined at Vt> — 0, which is the curve 
itself, and is numerically unstable near the curve. For this reason, we regard v 
as a distance function to a "tube" of small radius e around the curve, instead of 
extracting the true ID curve. That is, we evolve the e-level set instead of evolving 
the true curve (Figure 4(b)). Note that e does not denote a fixed value here: we 
mean simply that the evolving shape is a "tubular" surface of some (unspecified 
and variable) nonzero width. In addition to being more robust, this method 
better captures the geometry of blood vessels, which have nonzero diameter. 

Banding: Instead of evolving the entire volume, we evolve only the portion 
of the volume within a narrow band of the zero level set (the current surface). 
This technique is commonly used in level set methods. Normally, we set the 
band to include voxels that are up to 6 voxels away from the surface. We have 
increased this distance up to 12 for some experiments. The advantage of this 
technique is efficiency, and the disadvantage is that we may miss structures that 
are outside the band if the potential function g does not have a large enough 
capture range to attract the segmentation to these structures. This issue can be 
addressed by ensuring that g is compatible with the band size. 

Curvature Instead of Eigenvalues: For computational efficiency and be- 
cause of numerical instability of the gradient computations and thus the evolu- 
tion equation near Vw = 0, we remark that the level sets of the function v flow 
in the direction of the normal with velocity equal to the sum of their smaller 
principal curvature and the dot product of Vv with the image-based vector field 
d. Therefore, we compute the smaller curvature directly from v instead of as an 
eigenvalue of Psj v V 2 vPs/ v . 

Image Scaling: To control the trade-off between fitting the surface to the 
image data and enforcing the smoothness constraint on the surface, we add an 
image scaling term imscale to Equation 5 to obtain 
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v t = X(Vv(x, t), V v(x, t)) + imscale * — S7v(x, t) ■ H 



VI 
9 \W 

imscale is set by the user or can be pre-set to a default value. 



(6) 
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Gradient Directionality: Because vessels appear brighter than the back- 
ground, we weight the image term by the cosine of the angle between the normal 
to the surface and the gradient in the image. This cosine is given by the dot 
product of the respective gradients of v and I, so the update equation becomes 

v t = X(Vv(x, i), VVx, £)) + imscale * (Vv ■ VI) * — Vv(x, t) ■ H , — , (7) 

9 I v/ I 

For example, if the two vectors point in the same direction, then the brighter 
region is inside the surface and the darker region is outside; the angle between 
the vectors is 0, whose cosine is 1, so the image term is fully counted. However, 
if they point in opposite directions, the negative weighting prevents the evolving 
vessel walls from being attracted to image gradients that point in the opposite 
direction. 

Reinitializing Volume: As customary in level set segmentation methods, 
the volume v is periodically reinitialized to be a distance function: the zero level 
set S is extracted, then each point in the volume is set to be its distance to S. For 
our implementation, this reinitialization is itself a level set method. To obtain 
the positive distances, the surface is propagated outward at constant speed of 1, 
and the distance at each point is determined to be the time at which the surface 
crossed that point. A second step propagates the surface inward to obtain the 
negative distances analogously. For some experiments, we have used the Fast 
Marching Method [12] to implement these steps. 

Initial Surface: Figure 7 shows additional detail on the generation of the 
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Fig. 7. More detailed illustration of initialization part of algorithm. 



initial surface. This initial surface (and thus the initial volume) is normally 
generated by thresholding the MRA dataset. However, the method does not 
require that the initial surface be near the target surface but may use any initial 
surface. Figure 8 illustrates a vertical bar evolving into the segmentation of the 
first dataset in Figure 10. 
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Fig. 8. Illustration of a vertical bar evolving in a segmentation of the first dataset in 
Figure 10. 



Smoothing: As shown in Figure 7, the datasets may be pre-processed to 
reduce noise. For the results presented here, the raw datasets were convolved 
with an isotropic Gaussian of a = 0.5. 

Cleaning: We post-process the segmentations to remove any surface patches 
whose surface area is less than some threshold (a parameter of the method) to 
eliminate patches corresponding to noise in the original data. 

Vessel Radii Estimation: The larger principal curvature can be useful in 
measuring the radii of the vessels for a particular application, since radius is the 
inverse of curvature. This curvature can be easily computed when the smaller 
principal curvature is computed for the segmentation. We have added the option 
to color-code our segmentations based on vessel radii, as estimated from the local 
larger principal curvature of the tubular surface. 



7 Results 

We demonstrate segmentation results on four datasets, courtesy of the Surgi- 
cal Planning Laboratory, Brigham and Womens Hospital and Harvard Med- 
ical School (Figures 9, 10, and 11). All datasets had an initial resolution of 
.9375x. 9375x1. 5mm 3 (256x256x60 voxels). The final example only was resam- 
pled to .9375x.9375x. 9375mm 3 (256x256x96 voxels) before segmentation; the 
other segmentations were performed directly on the raw data. The images are 
not square (256x256) because uninteresting portions were cropped for efficiency. 
In Figure 10, the initial surface for the segmentation was a surface obtained by 
thresholding the raw dataset whereas in Figure 9 it was a tube as in Figure 8; 
imscale also varied as discussed below. For comparison, Figure 10 first shows 
results obtained by thresholding alone. Figure 11 shows an enlargement of a 
portion of the segmentations and corresponding maximum intensity projection 
considered in Figure 10. 

The following parameters were used in these experiments; all settings were 
chosen empirically. For our method, imscale varied across the datasets depend- 
ing on the noise present. A threshold Unit was used in Figure 10 to obtain the 
initial surface from the dataset; such a threshold was obviously not needed in 
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Fig. 9. The first image in each row is the maximum intensity projection of the raw data, 
and the second and third are the segmentation result from two orthogonal viewpoints. 
These results are obtained by our method where the initial surface was a vertical bar 
as showed in Figure 8. 
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Fig. 10. Results on three datasets are shown. For each image pair, the first image is the 
maximum intensity projection of the raw data, the second is the segmentation result 
from thresholding only and the third is the segmentation result using our method. 
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Figure 9. A cleaning threshold c indicated the minimum surface area of connected 
components of the surface to be retained in the post-processing "cleaning" step. 

For thresholding only, the threshold t t hresh was chosen and also the cleaning 
threshold c. For all datasets, ti n u was slightly higher than tthresh for the same 
dataset: although using a lower tthresh alone looks better after the cleaning step, 
the noise before cleaning worsened our results and led us to use a slightly higher 
value for initialization. 

Recall that obtaining the very small vessels is the goal of this work since the 
large vessels are easily segmented by thresholding. For this reason, imscale was 
set fairly high in the experiments in Figure 10 to obtain the small vessels, at 
the expense of also obtaining many imaging artifacts. A coarser segmentation is 
obtained in Figure 9 by choosing lower values for imscale. Although the results in 
this figure are only similar to those obtained by simple thresholding, the objective 
of the demonstration is academic: it shows that we capture the vasculature shape 
even when the initial guess is meaningless. In practice, better results are obtained 
using thresholding for initialization. 

When considering that the imscale parameter controls the trade-off between 
noise and small vessels in our method, and when comparing our method to 
thresholding alone, it is important to note that it would not be possible to 
similarly lower tthresh to obtain the small vessels (and noise) by thresholding 
alone. Lowering the threshold obtains large blobs in the volume which do not 
correspond to vessels. Our method is thus more powerful than thresholding alone. 

Finally, we demonstrate the capability to color-code the vasculature surface 
based on local curvature. Notice (Figure 12) that for a ribbon-like vessel, the 
flatter sides shows a large radius, and the sharply curved edges show a small 
radius. In this example, the colorscale is continuous from darkest to lightest 
intensities, with darkest indicating a radius of curvature < 1mm and lightest 
indicating a radius of curvature > 2mm. The curvatures output by our evolution 
have been smoothed by a 3x3x3 filter prior to coloring the surface. 

8 Future Work 

Vessels tend to appear thinner in our segmentations than in those obtained by 
thresholding. One possible reason is that our method uses gradients instead of 
intensities, so the vessel wall is found attracted to the strongest gradients, which 
may be fully inside the bright region indicated by thresholding. A second option 
is that the underlying mathematics of our algorithm assume that the vessels 
are ID curves, not tubular surfaces. We believe that our e-level set method 
allows the method to successfully handle tubular surfaces, but have not yet 
verified this analytically. A final potential reason for the discrepancy is that 
the segmentations obtained by thresholding may be thicker than the true blood 
vessels due to noise around the vessels. Future work will involve comparisons to 
manual segmentations which will provide ground truth to evaluate both methods. 
We also observe a lot of noise in our segmentations of the first and second 
datasets. As mentioned above, we could obtain much less noise at the expense of 
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Fig. 11. Enlargement of a portion of the second example from Figure 10. As above, 
the second image is the segmentation obtained by thresholding alone, and the third 
image is the result of our method. 




Fig. 12. Our method naturally allows estimation of local radii of curvature of the seg- 
mented vessels. In this image of a partial segmentation of the first dataset in Figure 10, 
the colorscale is continuous from darkest to lightest intensities, with darkest indicating 
a radius of curvature < 1mm and lightest indicating a radius of curvature > 2mm. 
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the thinnest vessels by lowering imscale. For the large amounts of noise in these 
datasets, noise is often indistinguishable from small vessels when only a small 
local neighborhood is considered, as in our algorithm. To address this problem, 
one could reduce noise prior to segmentation by filtering or incorporate a more 
sophisticated image measure into Equation 5. 

On the positive side, the segmentation of small vessels that were not ob- 
tainable by thresholding encourages us to continue in the development of this 
algorithm. Although still in preliminary stages, we believe that it has the poten- 
tial to yield effective segmentations of very thin vessels. 
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